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°f  the  transfer fu"ctl0"i  ^ 

Chapter  I of  this  report  is  devoted  to  a detailed  discussion  of  the  oblect 
res.oration  problem.  A solution  is  effected  via  the  method  of  singular 
value  decomposition  and  a representative  problem  is  solved  and  compared 
ith  a previous  solution  via  Tichonov  regularization.  Two  situations  ere 

up  ta OcSr  imdlSCUS8ad  “d  the  EI"al  ''e”1“s 


SUMMARY 


The  objective  of  this  program  is  thf  analysis  of  the  funda- 
mental limitations  placed  upon  object  restoration  by  the  presence 
of  noise  and  unavoidable  atmospheric  and  optical  system  degrada- 
tions. Our  approach  is  to  treat  individually  and  in  concert  the 
several  factors  which  limit  image  formation  (direct  problem)  and 
further  limit  image  restoration  (inverse  problem),  to  understand 
their  interrelation  and  order  of  importance.  Within  this  overall 
framework,  we  have  treated  various  problem  areas  comprising  the 
•chapters  of  this  report:  the  object  restoration  problem  viewed 

as  an  improperly  posed  inverse  problem,  sums  of  random  variables 
individually  governed  by  lognormal  probability  density  functions, 
and  moment  behavior  of  the  transfer  function  random  process,  in- 
cluding temporal  effects  and  those  due  to  random  amplitude  errors 
acting  in  conjunction  with  random  phase  errors,  realizability 
conditions  on  the  covariance  of  the  wavefront  aberration  function. 

hapter  I is  devoted  to  a detailed  discussion  of  the  object 
restoration  problem.  The  purpose  of  this  chapter  is  two-fold. 
First  to  discuss  the  peculiar  mathematical  behavior  of  the  object 
restoration  problem;  second  to  discuss  computational  techniques 
designed  to  handle  such  problems.  In  particular  the  method  of 
singular  value  decomposition  is  employed  to  effect  a numerically 
stable  inversion.  The  solution  obtained  via  singular  value 


iii. 


decomposition  is  compared  to  the  solution  obtained  previously  by 
Barakat  and  Blackman  via  Tichonov  regularization  when  the  image 
is  a bar  target  and  che  optical  system  is  an  aberration-free  slit 
aperture.  In  addition  to  the  image  being  noisy,  calculations  a^e 
also  displayed  for  the  important  case  where  both  image  and  point 
spread  function  are  noisy  (5%  noise  in  both). 

Research  conducted  during  the  first  part  of  the  funding 
period  was  basically  completed  at  that  time,  and  the  work  was 
written  up  in  final  form,  we  refer  the  reader  to  report:  RADC- 

TR-74-277 , October  197^  for  full  details.  The  research  was 
divided  into  three  topics: 

1.  Sums  of  independent  lognormally  distributed  random 
variables . 

2.  Realizability  conditions  on  the  covariance  of  the  wave- 
front  aberration  function. 

3 • Statistics  of  the  transfer  function:  temporal  and 

random  amplitude  effects. 

The  lognormal  probability  density  permeates  much  of  the  current 
literature  on  optical  propagation  through  the  turbulent  atmos- 
phere and  has  been  the  source  of  a good  deal  of  controversy  in 
various  contexts  of  the  general  problem.  Sums  of  lognormal  ran- 
dom variables  arise  through  aperture  averaging  and  other  related 


iv 


processes  . 


Because  the  lognormal  density  assumes  central  impor- 
tance in  many  atmospheric  propagation  problems,  certain  aspects 
of  the  behavior  of  (sums  of)  lognormally  distributed  variables 
were  considered.  Particular  results  include  the  derivation  of  an 
expression  for  the  characteristic  function  of  the  lognormal  prob- 
ability density,  which  is  then  used  to  calculate  densities  for 
sums  of  lognormally  districted  random  variables;  a simplified 
proof  that  the  lognormal  probability  density  is  not  determined  by 
ito  moments  even  thougn  they  exist;  an  investigation  of  the  ob- 
served "permanence"  of  the  lognormal  density  in  terms  of  its 
asymptotic  behavior,  leading  to  a soundly  based  explanation  of 
the  phenomenon.  Realizability  conditions  on  the  covariance  of 
the  wavefront  aberration  function  were  also  briefly  considered; 
the  class  of  admissible  covariance  functions  is  restricted  to 
those  representing  (isotropic)  random  surfaces.  Investigations 
of  the  statistics  of  tne  transfer  function  random  process,  pre- 
dominantly restricted  in  our  prior  work  to  errors  of  phase,  have 
been  e. [tended  to  include  the  statistical  dependence  of  che  trans- 
fer function  process  on  time  (and  concurrently  on  integration 
time),  as  well  as  on  jointly  present  random  amplitude  and  phase 
errors.  Results  show  that  the  expected  value  of  the  time  depen- 
dent transfer  function  does  not  depend  on  integration  time,  but 
that  the  general  higher  moment  behavior  is  collection  time 
dependent.  Random  amplitude  effects  are  considered  within  the 


usual  lognormal  and  normal  assumptions  regarding  the  respective 
amplitude  and  phase  probability  densities,  but  without  the  added 
restriction  of  independence  between  the  amplitude  and  phase  ran- 
dom processes.  In  terms  of  the  expected  transfer  function,  a 
product  of  four  terms  arises:  one  standard  deterministic  p^rt, 

one  part  due  to  phase  errors  acting  alone,  one  part  due  solely 
to  amplitude  perturbations,  and  a cross  term  which  is  of  unit 
modulus  and  depends  for  its  existence  on  the  odd  part  of  the 
cross  correlation  between  amplitude  and  phase. 
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CHAPTER  I 

THE  OBJECT  RESTORATION  PROBLEM  AS  AN 
IMPROPERLY  POSED  INVERSE  PROBLEM: 

A SOLUTION  VIA  SINGULAR  VALUE  DECOMPOSITION 
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1.  INTRODUCTION 

In  producing  an  image  of  a remote  object  through  an  inter- 
vening turbulent  atmosphere,  we  are  faced  with  the  reality  that 
the  image  is  an  imperfect  and  often  inadequate  representation  of 
the  object.  In  order  to  obtain  more  information  one  natural  ap- 
proach is  post-detection  compensation,  whereby  one  attempts  to 
obtain  the  object  given  the  image  and  the  relevant  diffraction 
characteristics  (i.e.,  point  spread  function  or  transfer  function) 
of  the  optical  system.  We  term  this  the  object  reconstruction 
problem.  Unfortunately  the  object  reconstruction  problem  is  one 
of  a class  of  inverse  problems  that  are  extremely  sensitive  to 
noisy  data;  in  fact  so  much  so  that  naive  and  ad  hoc  techniques 
are  virtually  powerless  to  effect  a viable  solution. 

This  fact  has  been  Implicitly  recognized  by  many  people  and 
the  existence  of  a current  program  on  real  time  predetecticn 
compensation  is  proof  of  this  realization.  Nevertheless,  even 
if  post-detection  compensation  is  currently  in  disfavor,  the 
fault  has  been  more  with  the  proporents  of  this  approach  rather 
than  with  the  approach  itself,  in  that  many  methods  currently 
advocated  and  implemented  on  a computer  are  not  capable  of 
handling  the  inherently  ill-posed  nature  of  the  problem. 

The  purpose  of  this  part  of  the  report  is  twc-fold.  First 
to  discuss  in  the  peculiar  mathematical  behavior  of  the  object 
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reconstruction  problem;  second,  to  discuss  computational  tech- 
niques speciiically  designed  to  handle  such  inverse  problems. 
Thus  this  part  of  the  report  contains  didactic  material  as  well 
as  original  material. 


It  is  probably  safe  to  say  that  inverse  problems,  such  as 
the  object  reconstruction  problem,  are  among  the  most  difficult 
problems  facing  the  scientific  community  today.  The  problems 
are  difficult  enough  in  themselves  and  the  overly  optimistic  at- 
titude of  the  optical  reconnaissance  community  that  brute  force 
computation  on  large  memory,  high  speed  computers  can  circumvent 
the  inherent  difficulties  has  tended  to  retard  progress  in  this 
complicated  subject.  So  much  for  the  polemics,  now  to  the  tech- 
nical problems. 


In  order  to  proceed  we  must  postulate  or  derive  a relation- 
ship becween  object  and  image.  Unless  we  possess  such  a mathe- 
matical or  physical  model  we  cannot  proceed.  If  the  optical  sys- 
tem is  illuminated  by  spatially  incoherent,  quasi-monochromatic 
illumination,  the  transfer  of  radiation  from  the  object  plane  to 
the  image  plane  is  a lirear  operation  in  the  illuminance.  Thus 
by  virtue  of  the  superposition  integral 


h (x , y ) 


00 

t(x,xf ;y,y ' )o(x ' ,y ' )dx'dy' 

— oo 


(1.1) 
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where : 


h(x,y)  = distribution  of  illuminance  in  the  image 
t(x>y)  = point  spread  function 
°(x,y)  = object  intensity  function. 

Equation  1.1  is  an  integral  equation  of  the  first  kind  for  the 
unknown  object  o in  terms  of  the  measured  (and  thus  noisy)  h and 

known  (but  possibly  noisy)  t.  We  consider  only  objects  )f  finite 
extent,  i.e., 


°(x,y)  = 0 for  x , y t R 


(1.2) 


where  R is  a finite  region  not  necessarily  symmetric  about  the 
origin  of  coordinates.  All  that  we  know  about  o is  that  it  lies 
in  the  larger  finite  region  A (i.e.,  o(x,y)  c:  R c A ) where  A is 
a region  about  the  origin  of  coordinates  for  which  we  have  mea- 
surements of  the  image.  A is  not  necessarily  symmetric,  but 
probably  would  be  in  most  application.  Thus,  Eq . 1.1  becomes 


:x,y)  = t(x,x';y,y»)o(x',y»)dx'dy'  . 


Note  that  we  are  not  postulating  the  isoplanatic  conditi( 

t(x,x';y,y')  = t (x-x ' ,y-y ' ) . 


(1.3) 


(1.4) 


There  are  no  restrictions,  such  as  symmetry,  placed  on  the  kernel 
t(x,y)  except  that  it  be  the  point  spread  function  of  an  optical 


system. 


Our  problem  Is  to  determine  the  unknown  o(x,y)  in  terms  of 
the  known  h(x,y)  and  t(x,y),  if  possible. 

-here  are  actually  two  topics  that  must  be  considered  in  an 
analysis  of  the  measurements  of  the  type  described  by  Eq . 1.1 
whicn  we  now  write  in  abstract  operator  form  as 

h = to  . (1.5) 

The  first  problem  is  the  mathematical  one  of  finding  an  inverse 
operator  t such  that  t 1 t is  the  unit  operator.  For  the  pur- 
pose of  the  present  discussion  such  an  inverse  is  assumed  to 
exist.  A "solution"  of  Eq . 1.5  is  then 

o = t~'h  . (1.6) 

The  solution  is  not  unique.  The  solution  o in  this  equation  is 
defined  only  to  within  an  additive  component  Oj  where  o)  is  any 
solution  of  the  homogeneous  equation 

to,  =0  • (1.7) 

A unique  solution  of  Eq . 1.1  is  possible  only  if  information 
sufficient  for  the  determination  of  o is  available.  The  mean- 
ing of  Eq.  1.7  is  that  the  measuring  system  (i.e.,  optical  sys- 
tem) is  completely  insensitive  to  those  spectral  components  of 
the  object  contained  in  Oj . However  since  we  are  going  to  deal 
only  with  finite  size  objects,  Eq . 1.2,  then  o = 0, 
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The  second  problem  in  the  analysis  of  Eq.  1.3  arises  from 
the  fact  tha^  both  h and  t are  subject  to  statistical  uncertain- 
ties. In  fact  small  errors  in  the  measured  image  h can  (and 
usually  are)  amplified  by  a nearly  singular  operator  t-1  to  such 
an  extent  that  the  "solution"  of  Eq . 1.3  is  completely  meaningless. 

The  second  problem  arises  from  the  fact  that  the  image  is  a 
measured  quantity  and  therefore  subject  to  statistical  uncertain- 
ties. It  will  turn  out  that  these  statistical  errors  are  extremely 
important.  Primary  causes  of  these  image  errors  are:  (1)  deter- 

ministic degradations  - aperture  and  aberration  limits;  (2)  sto- 
chastic degradations  such  as  residual  random  wavefront  errors  due 
to  the  combined  effects  of  turbulence  and  pre-detection  compensa- 
tion, Poisson  (at  best)  nature  of  return  and  Poisson  (or  similar) 
detector  noise.  It  is  tempting  to  consider  the  errors  as  arising 
from  additive  noise  so  that  h(x,y)  becomes 

hj (x,y)  + n(x,y ) (1.8) 

where  n(x,y)  is  the  space  variant  noise  and  hj(x,y)  the  noiseless 
image,  but  such  as  decomposition  is  simply  not  possible.  The 
reason  is  that  multiplicative  noise  also  enters  and  to  attempt  to 
separate  it  out  becomes  almost  Impossible.  It  is  better  to  con- 
sider the  image  h(x,y)  as  subject  to  noise  and  not  attempt  to 
specify  the  noise.  Also  the  point  spread  function  t(x,y)  is 
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generally  not  known  very  accurately.  Further  t(x,y),  especially 
for  the  types  of  problems  with  which  we  are  concerned,  can  even 
be  a random  function;  see  Barakat  (1971),  Barakat  and  Blackman 
(1973),  Barakat  and  Blackman  (197^). 

Given  this  background  information  let  us  see  why  the  inverse 
operator  t-1  is  so  badly  behaved.  The  mathematical  relation  be- 
tween h and  o given  by  Eq.  1.3  (or  equivalently  Eq . 1.5)  is  basic- 
ally a smoothing  operation  if  we  view  o as  given  and  seek  to  cal- 
culate h.  If  the  optical  instrument  were  perfect,  then  assuming 
the  isoplanatic  condition,  we  would  have 

t(x,x* ; .V , y ' ) - 6 (x-x ’ ) 6 (y-y ' ) . (1.9) 

Thus  Eq.  1.3  would  become 


h (x , y ) = o(x,y ) , x,y  e A 

= 0 , x , y t A . (1.11) 

Of  course,  any  point  spread  function  satisfying  Eq.  1.10  is  not 
physically  realizable.  At  best  one  must  settle  for  a broadened 
response,  which  has  a maximum  in  the  vicinity  of  x = x',  y = y' 
and  tails  off  to  zero  as  x - x*,  y-y'  -*■  ±°°.  Such  a point  spread 
function  necessarily  smooths  the  object  somewhat  causing  a loss 


7 


f 


°f  information  in  the  measured  image  h.  This  can  be  seen,  for 
example,  more  normally  via  the  Riemann-Lebesque  theorem,  which 
states  that  for  irtegrable  t 
lim  r (• 

a 8 * 00  JJ  t(x,x';y,y')  elax ' e13y'  dx'dy*  = 0 . (1.12) 

A 

xhus  an  arbitrarily  high  frequency  component  of  o(x,y)  has  a small 
effect  on  h(x,y).  This  is  just  another  way  of  saying  that  the 
optical  system  has  a bandlimited  point  spread  function. 

Our  problem  is  to  recover  this  information  and  thereby  deter- 
mine the  unknown  object  o(x,y);  we  term  this  inversion  problem  the 
object  restoration  problem.  This  is  a classical  ill-posed  prob- 
lem in  the  sense  of  Hadamard . According  to  Hadamard,  a problem 
is  well  posed  if  the  following  three  conditions  are  satisfied: 

(a)  solution  exists 

(b)  solution  is  unique 

(c)  solution  depends  continuously  on  the  input  data. 

We  need  not  concern  ourselves  with  conditions  a and  b as  they 
are  satisfied  in  the  object  restoration  problem.  Barakat  (un- 
published) has  given  formal  proofs  of  conditions  a and  b.  How- 
ever condition  c is  not  satisfied  as  evidenced  by  the  Riemann- 
Lebesque  theorem,  Eq.  1.1?;  object  restoration  problems  do  not 
satisfy  condition  c since  a small  change  in  the  image  data  h can 
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correspond  to  an  arbitrarily  large  change  in  the  solution  c. 

Thus  the  inverse  operator  f1  does  not  have  a bounded  inverse. 
This  numerical  instability  is  inherent  in  the  very  nature  of  the 
problem. 

Our  basic  problem  is  to  stabilize  the  solution  against  the 
numerical  instabilities  of  such  ill-posed,  problems. 
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2.  ERROR  ANALYSIS  AND  CONDITION  NUMBER 

Given  the  preliminary  information  of  the  previous  sections, 
we  now  undertake  to  discuss  in  some  detail  the  technical  problems 
which  must  be  faced  in  attempting  object  reconstruction  on  a 
rational  basis. 

To  begin  with,  we  might  as  well  face  the  fact  that  the  con- 
tinuous version  of  the  problem  is  basically  useless  since  the 
image  is  measured  at  discrete  points.  Therefore  we  must  replace 
the  continuous  version  by  a discrete  version.  The  double  integral 
is  replaced  by  a quadrature  formula  (the  actual  one  employed  is  of 
no  great  concern)  and  the  continuous  variables  by  discrete  mesh 
points.  If  the  weights  for  the  quadrature  formula  are  denoted 
by  Hk  and  the  quadrature  points  by  xR,  etc.,  we  obtain 


n n 


h(xi-yj)  ' > ;2-1> 

where  i,j=l,***,m.  If  we  employ  lexicographic  ordering  we  can 
write  this  more  concisely  In  matrix  notation  as 


A A A 

h = to  , 


(2.2) 


The  matrix  t is  m x n (m  rows,  n columns)  and  in  general  m > n; 

/V  A 

h and  o are  in  this  concise  rotation  column  vectors  of  size  m 
and  n respectively. 
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The  validity  of  this  discretization  depends  on  two  approxi- 
mations. First,  the  replacement  of  an  infinite  dimensional  func- 
tion space  by  a finite  dimensional  vector  space.  Second,  the 
replacement  of  an  integral  operator  by  a finite  matrix.  It  will 
be  assumed  that  both  approximations  are  adequate  and  attention 
will  be  denoted  to  solving  the  resultant  system  of  linear  equations, 
Eq.  2.2. 

The  matrices  resulting  from  ill-posed  problems,  such  as 
Eq.  1.3j  are  inherently  ill-conditioned  irrespective  of  the 
particular  discretization  (quadrature)  scheme  employed,  so  long  as 
the  quadrature  scheme  is  reasonably  faithful  in  approximating  the 
integral . 

Let  us  consider  the  sensitivity  of  the  solution  of  our  canoni- 

A A 

cal  set  of  equations,  Eq.  2.2,  to  small  variations  in  t and  o 
where  m = n so  that  t is  square.  We  write 

t(o+6o)  = fi  + 6fi  , (2.3) 

where  6o,  6h  are  perturbations  of  the  object  and  image.  Then 

t6o  = 6h  , ( 2*  *0 

or 

6o  = t_ 1 6h  . (2.5) 
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It  is  convenient  to  have  a single  number  which  gives  an  overall 
assessment  of  the  "size"  of  a matrix  and  plays  the  same  role  as 
the  modulus  in  the  case  of  a complex  number.  For  this  purpose 
the  Euclidean  norm  is  introduced 

ll?ll  = <Z  i tjj’)*  . (2.6) 

By  Schwartz's  inequality 

I 1 6° ||  = lit"1  6h|  | < | It"1  | | | |6hj  | , (2.7) 

we  have 


jiiiii  < ||t- 

| 1 6h  | 


(2.8) 


For  the  relative  change  in  the  image,  | | 5o | |/ | | 6 | | ; we  take  the 

norm  o ^ Eq . 2.2  and  employ  Schwartz's  inequality,  with  the  final 
result  being 


o > h 


(2.9) 


Combining  unis  result  with  Eq.  2.3,  yields 

lliiil  < llt-'ii  'Mfili  - m£|| 

1151!  “ I |t| |->  ||h|| 


This  equation  shows  that  the  relative  change  in  the  object  due  to 
relative  change  in  the  observed  image  depends  upon  the  quantity 
l|t||  ||t-1||.  We  will  term  this  quantity  a condition  number , 

.A 

k(t),  to  indicate  how  the  change  of  the  object  solution  depends  on 


the  change  in  the  input. 


Since  we  have  already  indicated  that  tht  object  restoration 
problem  is  improperly  posed,  we  must  expect  that  the  condition 
number  k(t)  will  be  very  large.  For  example  if  k(t)  = 106  (a  com' 
mon  value  in  object  restoration  problems),  a perturbation  of  2“ 12 
in  the  elements  of  t can  change  the  computed  solution 

0 = h , (2.11) 

by  a factor  of  106  2" 1 2 = (5/2) 6 , that  is  even  the  leading  digit 
is  wildly  inaccurate!  A theoretical  upper  bound  for  the  rexative 
error  in  the  computed  solution  is  given  by 


I I 6o I 1 < k(t) 

ll°ll  “ l-k(t)J-bi 

lit 

provided 


lii?.  I 1 + Jiiili 
I Ih|  | ||t|| 


(2.12) 


A 

| |6t  | | 


(2.13) 


This  theorem  gives  the  upper  bound  of  the  '’■ariation  of  o due  to 
perturbation  in  t and  h.  Thus  ill-conditioning  is  associated  with 
a large  value  for  k(t),  since  in  this  case  the  bound  on  the  error 

A 

in  o is  very  large.  This  behavior  corresponds  to  the  unboundedness 
of  the  inverse  of  the  integral  operator  from  which  t is  obtained. 


Thio  discouraging  state  of  affairs  is  characteristic  of 
inversion  problems  and  we  must  accept  the  fact  that  naive  attempts 
at  the  solution  of  the  problem  such  as  direct  inversion,  Eq.  2.11 


13 


will  not  succeed.  Unfortunately  this  fundamental  fact  has  not 
been  realized  by  most  previous  investigators  who  have  placed 
undue  reliance  on  machine  computation  under  the  assumption  that 
a larger  matrix  t will  lead  to  better  results. 

Since  direct  inversion  will  not  work,  it  is  tempting  to 
consider  least  squares  as  an  inversion  technique. 


3.  LEAST  SQUARES  INVERSION 

If  m>n,  then  the  natural  method  of  solution  that  parallels 
direct  inversion  (when  the  matrix  t is  square)  is  to  minimize 

|to  - hj I2  (3.1) 

This  is  equivalent  to  solving  the  linear  system 

A A 4>  A A 

t h - t to  (3.2) 

A 

where  t denotes  the  transpose  of  t.  The  matrix  of  coefficients 
t+t  is  now  symmetric.  The  o which  satisfies  Eq.  3.2  is  called 
the  least  squares  solution  and  is  given  by 


A 

o 


(t+er>t+h 


(3.3) 


However,  in  using  Eq.  3.3  it  is  very  possible  that  t+t  may  be 

A 

even  more  ill-conditioned  than  t itself!  Thus  an  attempt  to 

A | A 

invert  the  square  matrix  (t  t)  will  produce  meaningless  results 
similar  to  those  obtained  by  direct  inversion.  The  reason  for 
this  paradoxical  state  of  affairs  is  bound  up  with  the  fact  that 

A.J.  A 

t t is  usually  rank  deficient  even  though  it  is  formally  over- 
determined in  that  there  are  more  columns  than  rows.  That  is  to 
say,  the  rank  should  be  determined  during  the  course  of  computing; 
unfortunately  information  about  rank  deficiency  cannot  be  obtained 
from  triangular  factorization  as  employed  in  least  square 
calculations . 
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Mdiniia 


m 


Now  one  way  to  build  up  the  rank  deficiency  is  to  augment 
the  data  provided  by  the  optical  system  (u,  the  image)  with  addi- 
tional a priori  knowledge  of  the  nature  of  the  physical  problem 
in  order  to  make  the  computed  solution  at  least  physically  mean- 
ingful. Such  constraints  usually  take  the  form  of  inequalities 
imposed  on  the  solution.  These  constraints  can  be  either  implicit 
or  explicit.  The  task  is  to  select  from  the  infinite  number  of 
possible  solutions  which  satisfy  the  observational  data  within 
experimental  error,  the  one  which  best  satisfies  some  set  of 
implicit  or  explicit  constraints.  We  begin  by  discussing  the 
Tichonov  regularization  algorithm  which  employs  implicit  constraints 
and  which  was  employed  by  Barakat  and  Blackman  (1973)  in  their 
studies  on  object  restoration. 
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4.  TICHONOV  REGULARIZATION  ALGORITHM 

Most  natural  objects  are  smooth.  This  vague  intuitive  con- 
cept has  given  rise  to  a family  of  inversion  algorithms  which 
require  the  solution  to  be  smooth  in  some  prescribed  sense.  These 
algorithms  go  under  the  general  name  of  regularization  and  are 
associated  with  the  mathematician  Tichonov. 

The  basic  idea  behind  regularization  is  the  replacement  of 
the  least  squares  minimization,  Eq.  3-l>  by  minimization  of  the 
functional 

Ma  = ||tS-h||2  + a | | o ' | | 2 

where  a,  termed  the  regularization  parameter,  is  a small  non- 
negative number.  It  can  be  shown  that  the  solution  via  this 
approach  is  numerical  stable  in  the  norm  sense.  The  functional 
Ma  is  minimized  by  the  usual  machinery  of  the  calculus  of  varia- 
tions by  setting  the  first  variation  equal  to  zero.  We  omit  the 
details  since  they  are  available  in  Barakat  and  Blackman  (1973). 
The  final  result  is  that 

o = (t  t + aH)  1 t h (4.2) 

A 

where  A is  the  matrix  of  central  difference  operators.  If  a is 
large  [i.e.,  a ~ o(l)]5  then  the  regularization  term  tends  to 
swamp  the  first  term  and  the  object  solution  vector  is  much  too 
smooth.  On  the  other  hand,  if  a = o then  the  equation  reduces  to 


17 


1 


o 


A A A - A A A 

(t+t)_1  t+h 


( ^.3) 


which  is  numerically  unstable,  being  the  usual  least  squares 
solution  already  discussed.  The  choice  of  a depends  on  the  shape 
of  the  object  being  reconstructed  and  on  the  noise  level  in  the 
image . 

Barakat  and  Blackman  (1973)  have  utilized  the  Tichonov 
regularization  algorithm  in  their  object  reconstruction  studies. 
In  the  next  section  we  will  discuss  the  use  of  singular  value 
decomposition  for  object  reconstruction  and  in  order  to  compare 
both  methods,  we  first  summarize  some  of  the  Barakat-Blackman 
work.  Their  numerical  calculations  (not  the  theory)  is  confined 
to  that  of  a one-dimensional  unit  pulse  of  half  width  xQ  imaged 
by  an  aberration-free  slit  aperture.  Thus  the  object  and  the 
point  spread  function  are  given  respectively  by 


o(x)  = 0 , -°°  < x < — xQ 

= 1 , -X  < X < X 

’ o o 

= 0,  x < x < °°  ( *4 . *0 

’ 0 


and 


with  the  normalization  constant  tt-1  determined  by  the  constraint 
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t(x)dx  = 1 . 


(4.6) 


We  now  assume  that  the  isoplanatic  condition  holds  so  that 


Eq.  1.3  reads 


* 

l(x  ) = t (x-x  ' )o  (x  ' )( 
A 


where  A is  the  integral  (-10  < x < 10).  All  we  know  is  that 
o(x)  lies  within  this  interval  as  discussed  in  Sec.  1. 

To  compute  the  image  of  the  pulse,  we  note  that 


h (x ) = ± 


(**&■)  • 


x0  < 10 


This  can  be  evaluated  explicitly  in  terms  of  the  sine  integral 


Si(z)  = 


sinx 


so  that 


h(x)  = j jSi[2(x+xo)  - Si [2 (x-x  )] 


(x+x  ) sin2 (x-x  ) J 

TT ~ V + 


,X  + X0  )' 


The  image  h(x)  for  x = M is  shown  in  Figs.  1 and  2 (see  the 


f 


dotted  lines),  while  the  object  o(x)  given  by  Eq . 4 . 4 is  the 
heavy  solid  line. 

Barakat  and  Blackman  inverted  Eq . 4.7  via  the  Tichonov 
algorithm  for  two  cases:  one,  the  academic  noise  free  case  shown 
in  Fig.  1;  two,  the  noisy  case  (5$  noise  in  image)  shown  in  Fig. 
2.  The  necessary  analytical  and  computational  details  are  given 
in  their  paper  along  with  numerical  data  on  other  situations. 

Note  that  in  applying  the  Tichonov  algorithm,  we  have  made 
practically  no  use  of  any  a priori  information  about  the  object 
except  that  it  lies  somewhere  in  the  interval  10  < x < 10.  The 
reconstructed  object,  shown  in  the  various  figures,  has  the  un- 
realistic feature  of  negative  illuminance  over  small  intervals. 
However,  considering  that  we  did  not  demand  non-negativity  of 
the  reconstructed  ooject  the  results  are  very  good  especially  for 
the  noisy  case . 


20 


5. 


SINGULAR  VALUE  DECOMPOSITION 
Although  the  Tichonov  regularization  algorithm  Is  very  power- 
ful, the  fact  that  the  resulting  object  reconstructions  can  be- 
come negative  is  a mark  against  the  method.  Therefore  we  turn 
to  another  inversion  algorithm  which  is  even  more  powerful,  the 
method  of  singular  value  decomposition.  It  should  be  pointed  out 
that  the  principal  investigator  has  recently  utilized  this 
algorithm  to  invert  photoelectron  correlation  function  data  to 
obtain  spectral  line  shapes  (Barakat  and  Blake,  1975). 

The  singular  value  decomposition  of  the  m * n (m  £ n)  real 
matrix  t is  given  by  the  factorization 


A A A A 1 

t = UDV 


(5.1) 


where  U is  an  m * m orthogonal  matrix  and  V is  an  n * n ortho- 
gonal matrix 


AiA  A Al 

u u = UU  = I 


m 


A X A A Al  A 

V V = W = In  . 


(5.2) 


Here  D is  an  m * n matrix  whose  only  nonzero  elements  are  on  the 
principal  diagonal 


D » diag(o1,o2,***on,  0,0, •••,0) 


(5-3) 


where 
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°1  - a2  - 


(5.4) 


and  the  remaining  (m  - n)  diagonal  elements  are  zero. 


The  columns  of  U can  be  shown  to  be  the  orthonormal  eigen- 
vectors of  t+t,  while  the  columns  of  V are  the  orthonormal 

eigenvectors  of  tt+.  Finally  the  singular  values  of  t,  o^,  are 

J 

mathematically  equal  to  the  non-negative  square  roots  of  the 
symmetric  matrices  tt+,  t+t 


^ « A A A 

o = yltt  o 


A n A ^ A A 

o = yit  to 


(5.5) 


so  that 


0,  = +(v2.)^  . 

d d 


(5-6) 


However  a word  of  caution;  the  tact  that  Oj  and  y are  related 
in  such  a simple  manner  might  tempt  one  to  believe  that  the  prob- 
lem is  perfectly  straightforward  since  one  merely  solves  the 
eigenvalue  problem  stated  in  Eq.  5.5.  However,  singular  values 

A 

correct  to  working  accuracy  for  t can  often  be  computed  when 
certain  small  eigenvalues  cannot  be  computed  for  t+t  or  tt+.  To 
anyone  who  has  ever  done  serious  computing  this  fact  is  not 
startling;  it  is  caused  by  the  perturbation  of  an  exact  t+t 
introduced  in  the  multiplication  of  t+  by  t. 
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The  condition  number  k(t)  of  t can  be  expressed  in  terms  of 


its  singular  values 


k ( t ) = 


where  o and  o . are  the  maximum  and  minimum  singular  values 
of  t.  Thus  a ill-conditioned  matrix  is  one  with  a great  varia- 
tion in  the  magnitude  of  its  singular  values.  Values  of  k(t)  = 
o(106)  are  common  (i.e.,  have  been  encountered  by  the  author 
during  the  course  of  the  numerical  work).  To  see  how  this  af- 
fects the  "solution"  of  our  problem  h = to,  let  us  substitute 
Eq . 5.1  into  Eq.  2.2.  Upon  performing  the  calculations  we  have 


A A A+A  + A 

o = VD  U h 


where 


= dlag( cj“  1 , • • *0” 1 , 0 , • • • , 0 ) . 


The  expansion  becomes  clearer  if  the  summation  is  written  out 


explicitly 


u,h 

A n 1 A 

O = ) — ^ — V, 

J °J  J 


here  uj>vj  are  the  column  vectors  of  U and  V respectively. 

The  smaller  singular  values,  entering  into  the  denominator  tend 
to  magnify  greatly  any  error  in  the  measured  image  data  vector  h 
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resulting  in  a spurious  solution.  To  alleviate  this  state  of  af- 
fairs, we  must  cut  off  the  expansion  before  the  contamination  due 
to  the  small  singular  eigenvalues  creeps  in.  Specifically  let  us 
set 

0“*  = 0 if  °j  > e (5.11) 

where  a reasonable  criterion  for  picking  e is 

^— >>  noise.  (5.12) 

0 

In  applying  the  singular  value  decomposition  to  our  object  re- 
storation problem  we  employed  the  programming  techniques  described 
in  the  basic  paper  by  Golub  and  Reinsch  (1970). 

In  order  to  demonstrate  the  power  of  the  singular  value  de- 
composition algorithm,  we  again  consider  the  problem  discussed 
in  the  previous  section.  The  results  of  the  calculations  are 
summarized  in  Figs.  3 and  In  each  case,  the  solution  comprises  a 

data  set  on  the  interval  |x|  < 10.  Figure  3 is  to  be  compared 
with  Fig.  1,  and  Fig.  4 with  Fig.  2.  Two  facts  are  pertinent. 

One,  the  singular  value  reconstruction  has  only  very  small  nega- 
tive intensities;  two,  the  slopes  of  the  singular  value  recon- 
struction are  much  higher  than  those  of  the  Tichonov  algorithm 
even  for  the  noisy  case. 
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6. 


NOISY  POINT  SPREAD  FUNCTIONS 


I 


Thus  for  the  integral  operator  t (or  equivalently  the  matrix 

A 

t)  has  been  assumed  to  be  subject  to  no  uncertainties.  However, 
as  we  have  already  pointed  out,  the  point  spread  function  in  many 
situations  is  itself  subject  to  noise  and  it  is  important  to  assess 
the  influence  of  small  errors  on  t on  the  solution  o.  Ideally 
one  would  like  to  know  under  what  conditions  the  solution  o 
produced  by  a given  t is  a continuous  function  of  the  matrix 

A 

elements  of  t,  as  well  as  a quantative  measure  of  the  possible 
effects  of  matrix  errors. 

Some  aspects  of  this  difficult  problem  have  already  been 
investigated  by  Barakat  and  Blackman  (1973  )•  They  studied  the 
direct  problem 

h = ° (6.1) 

where  t is  a random  operator  for  the  edge  spread  function  and 
calculated  the  expected  value  of  the  edge  spread  function. 

In  view  of  this  fact,  it  was  felt  that  some  calculations 
should  be  made  of  the  Influence  of  a noisy  6 on  object  reconstruc- 
tion. The  author  made  an  attempt  to  quantify  -che  influence  of 
noisy  t on  o,  but  was  not  very  successful  since  the  bounds  on  o 
were  only  weakly  related  to  5t. 
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Therefore  it  was  decided  that  some  numerical  calculations 

A /V 

involving  both  a noisy  t and  h would  be  useful.  As  before  we 
confine  ourselves  to  the  situation  described  by  Eqs.  4.4  and  4.5> 
singular  value  decomposition  was  employed.  Both  the  point  spread 
function  and  the  image  wer^  subjected  to  5%  noise.  The  results 
of  a typical  object  reconstruction  are  illustrated  in  wig.  5* 

The  reconstructed  object  is  now  slightly  asymmetric  die  to  the 
fact  that  the  noisy  point  spread  function  is  no  longer  symmetric. 
Even  with  the  added  burden  of  a noisy  t,  it  would  appear  that  the 
reconstructed  object  via  singular  value  decomposition  is  "better" 
than  the  reconstructed  object  via  Tichonov  regularization;  compared 
Pigs.  5 and  2. 

Thus  is  would  appea  • that  singular  value  decomposition  offers 
a viable  approach  to  object  reconstruction  when  both  point  spread 
function  and  image  are  mildly  corrupted  by  noise.  The  fact  that 
the  reconstructed  object  has  a small  amount  of  negative  illuminance 
does  not  appear  to  be  a serious  problem. 
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FIGURE  LEGENDS 


Pig.  1.  Object  reconstruction  by  21  point  Tichonov  regularization: 
original  object,  — image,  • reconstructed  object. 

Pig.  2.  Object  reconstruction  by  21  point  Tichonov  regularization: 
original  object,  — — — image  (5#  noise  added), 

• reconstructed  object. 

Pi£*  3*  Object  reconstruction  by  singular  value  decomposition: 

heavy  solid  line,  object;  regular  solid  line,  reconstructed 
object;  dotted  line,  image. 

Fig.  Object  reconstruction  by  singular  value  decomposition; 

heavy  solid  line,  object  with  5%  noise  added;  regular 
solid  line,  reconstructed  object;  dotted  line,  image. 

Fig.  j.  Object  '^construction  by  singular  value  decomposition 

when  both  point  spread  function  and  image  are  corrupted 
by  5$  noise:  heavy  solid  line,  object;  regular  solid 

line,  reconstructed  object;  dotted  line,  image. 
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APPENOIX  A 

The  research  funded  during  the  beginnings  of  this  contract 
are  written  up  in  full  detail  in  our  previous  report:  RADC-TR- 

74-277,  October  1974.  Since  this  work  was  basically  completed  at 
that  time,  we  refer  the  reader  to  that  report  for  full  details. 

The  research  was  divided  into  three  topics,  they  are: 

1.  Sums  of  Independent  Lognormally  Distributed  Random 
Variables . 

2.  Realizability  Conditions  on  the  Covariance  of  the 
Wavefront  Aberration  Function. 

3.  Statistics  of  the  Transfer  Function:  Temporal  and 

Random  Amplitude  Effects. 
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The  following  errata  were  noted  in  the  interim  report: 

On  p.  3,  Eq.  (2),  the  denominator  should  contain  an  additional 


factor  of  2 inside  the  square  root. 

On  p.  last  paragraph,  for  x,,/x_  = y,  read  x,  /x  = ey  . 

K O K O 

On  p.  5,  Eq.  (10),  an  additional  factor,  elxy,  is  required  in 
the  integrand. 

On  p.  8,  top  line,  for  On  the. . . , read  One  of  the. . . 

On  p.  9,  Eq.  (24),  for  , read  (u  /u  2)  - 3 . 

2 **  2 

^2 

and  subtract  3 from  the  quotient  on  the  right-hand  side  of 

the  first  line;  the  second  line  is  correct  as  is 

On  p.  12,  first  of  Eqs.  (27),  denominator  should  contain  . 

On  p.  13,  Eq.  1,32) , the  factor  multiplying  t4  in  the  expansion 

of  the  logarithm  should  have  a numerator:  \i  - 3a4,  and 

<♦ 

the  remainder  is  0(N~3/2),  not  0(N_5/2). 

On  p.  13,  Eq.  (33),  for  (s  < 1),  read  (|s|  < 1)  . 

On  p.  1*4,  Eq.  (31*),  the  fourth  term  in  brackets  goes  as  t6,  not  t4. 

On  p . 1*4,  Eq . (37)>  the  second  term  in  brackets  should  have  a 

plus  sign  in  front. 

On  p.  15,  line  3>  for  tensive,  read  tensively. 
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On  p.  15,  second  paragraph,  line  8,  read  "the  degree  of 

approximation  is  governed  by  the  term  proportional  to  the 
coefficient  of  excess..."  (The  coefficient  of  excess  is,  by 
definition,  independent  of  N. ) 

On  p.  16,  Eqs.  (40)  and  (4l)  are  v.  f the  wrong  sign.  (See 
comments  above  on  Eq.  (37).)  The  sentence  beginning, 

"The  reason  for  the  negative  sign..."  shoald  therefore  be 
deleted.  It  is  Intuitively  obvious  that  the  skewness  of  a 
lognormal-like  distribution  must  be  positive.  (In  the 
penultimate  sentence  in  that  paragraph,  for  recultant, 
read  resultant . ) 

On  p.  18,  lines  8-9,  read  need  an  algorithm  which  retains.... 

On  pp.  23,  24,  the  standard  deviation  is  0.25,  not  the  variance. 

On.  p.  4l,  Eq.  (21),  the  z's  on  the  right  hand  side  of  equation 
should  be  lower  case. 

On  p.  43,  Eq.  (24),  in  the  third  term  on  the  right  hand  side  of 
the  equation  replace  by  a in  the  fourth  term  replace 


r . • by  r . , 
t^w  wiJj 


36 


